You can complete all the coding exercises in the JAGS Primer using simple R scripts. However, you might be wondering about how we created the JAGS Primer key you are looking at right now. We did this using Yihui Xie’s knitr package in R. This can be a highly useful tools for organizing your Bayesian analyses. Within the same document, you can:
Best of all, knitr can produce beautiful html files as output, which can be easily shared with collaborators. We encourage you to become familiar with knitr. We recommend Karl Broman’s knitr in a nutshell as an excellent introductory tutorial.
There is no x in the posterior distribution in equation 4. What are assuming if x is absent? Draw the Bayesian network, or DAG, for this model. Use the chain rule to fully factor the joint distribution into sensible parts then simplify by assuming that \(r\), \(K\) and \(\tau\) are independent.
There is no \(x\) because we are assuming it is measured without error.
Factoring the joint using the chain rule: \[[r,K,\tau, \mathbf{y}]=[y\mid r,K,\tau][r|K,\tau][K|\tau][\tau]\]
Simplifying under the assumption of \(r,K,\tau\) are independent: \[[r,K,\tau, \mathbf{y}]=[y\mid r,K,\tau][r][K][\tau]\]
We would obtain the same result by inspecting the DAG and using the rules:
A recurring theme in this course will be to use priors that are informative whenever possible. The gamma priors in equation 4 of the JAGS Primer include the entire number line \(>0\). Don’t we know more about population biology than that? Lets, say for now that we are modeling the population dynamics of a large mammal. How might you go about making the priors on population parameters more informative?
A great source for priors in biology are allometric scaling relationships that predict all kinds of biological quantities based on the mass organisms (Peters, 1983; Pennycuick,1992). If you know the approximate mass of the animal, you can compute broad but nonetheless informative priors on \(r\) and \(K\). This might leave out the social scientists, but I would trust the scaling of \(r\) for people if not \(K\).
In the absence of some sort of scholarly way to find priors, we can at least constrain them somewhat. There is no way that a large mammal can have an intrinsic rate of increase exceeding 1 – many values for \(r\) within gamma(.001, .001) are far large than than that and hence are complete nonsense. We know \(r\) must be positive and we can put a plausible bound on its upper limit. The only requirement for a vague prior is that its “\(\ldots\) range of uncertainty should be clearly wider that the range of reasonable values of the parameter\(\ldots\)” (Gelman and Hill, 2009, page 355), so we could use \(r\) ~ uniform(0, 2) and be sure that it would be minimally informative. Similarly, we could use experience and knowledge to put some reasonable bounds on \(K\) and even \(\sigma\), which we can use to calculate \(\tau\) as \(\tau=\frac{1}{\sigma^{2}}\).
Peters. The ecological implications of body size. Cambridge University Press, Cambridge, United Kingdom, 1983.
C. J. Pennycuick. Newton rules biology. Oxford University Press, Oxford United Kingdom, 1992.
A. Gelman and J. Hill. Data analysis using regression and multilievel / hierarchical modeling. Cambridge University Press, Cambridge, United Kingdom, 2009.
for loopsWrite a code fragment to set vague normal priors for 5 regression coefficients – dnorm(0, 10E-6) – stored in the vector b.
for(i in 1:5){
b[i] ~ dnorm(0, .000001)
}
Write R code (algorithm 3) to run the JAGS model (algorithm 3) and estimate the parameters, \(r\), \(K\) \(\sigma\), and \(\tau\). We suggest you insert the JAGS model into this R script using the sink() command as shown in algorithm 4. because this model is small. You will find this a convenient way to keep all your code in the same R script. For larger models, you will be happier using a separate file for the JAGS code. Here is the joint distribution for our logistic model again, with the priors updated from exercise 2 and \(\tau\) expressed as a derived quantity:
\[\begin{eqnarray} \mu_{i} & = & r-\frac{rx_{i}}{K}\textrm{,}\nonumber\\[1em] \tau & = & \frac{1}{\sigma^{2}}\textrm{,}\nonumber\\[1em] \left[r,K,\sigma\mid\mathbf{y}\right] & \propto & \prod_{i=1}^{n}\textrm{normal}\left(y_{i} \mid \mu_{i},\tau\right)\textrm{uniform}\left(K\mid 0,4000\right) \textrm{uniform}\left(\sigma\mid 0, 5\right) \textrm{uniform}\left(r\mid 0,2\right)\textrm{.}\nonumber\\ \end{eqnarray}\]
We use the sink() command to create a JAGS script from our joint distribution. This file is created within R and saved in the working directory. Please note that the outer set of brackets are only required when running this code within an R markdown document (as we did to make this answer key). If you are running them in a plain R script, they are not needed.
{ # Extra bracket needed only for R markdown files
sink("LogisticJAGS.R")
cat("
model {
# priors
K ~ dunif(0, 4000)
r ~ dunif (0, 2)
sigma ~ dunif(0, 50)
tau <- 1/sigma^2
# likelihood
for(i in 1:n) {
mu[i] <- r - r/K * x[i]
y[i] ~ dnorm(mu[i], tau)
}
}
",fill = TRUE)
sink()
} # Extra bracket needed only for R markdown files
Then we run the remaining commands discussed in the JAGS Primer. Note that jm is calling the JAGS script LogisticJAGS.R we just created.
#Be sure to do this sorting. It will be important for plotting later.
Logistic <- Logistic[order(Logistic$PopulationSize),]
inits = list(
list(K = 1500, r = .2, sigma = .01),
list(K = 1000, r = .15, sigma = .5),
list(K = 900, r = .3, sigma = .01))
data = list(
n = nrow(Logistic),
x = as.double(Logistic$PopulationSize),
y = as.double(Logistic$GrowthRate))
n.adapt = 1000
n.update = 10000
n.iter = 10000
jm = jags.model("LogisticJAGS.R", data = data, inits = inits, n.chains = length(inits), n.adapt = n.adapt)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 3
## Total graph size: 212
##
## Initializing model
update(jm, n.iter = n.update)
zm = coda.samples(jm, variable.names = c("K", "r", "sigma", "tau"), n.iter = n.iter, n.thin = 1)
MCMCsummary(zm)
## mean sd 2.5% 50% 97.5% Rhat
## K 1.237502e+03 6.270432e+01 1.128602e+03 1.232806e+03 1.375907e+03 1
## r 2.006954e-01 9.686011e-03 1.813155e-01 2.006520e-01 2.196546e-01 1
## sigma 2.864351e-02 3.040275e-03 2.345873e-02 2.839627e-02 3.531073e-02 1
## tau 1.259338e+03 2.612025e+02 8.020225e+02 1.240159e+03 1.817152e+03 1
## n.eff
## K 7016
## r 7235
## sigma 15059
## tau 16330
K, one as a smooth curve and the other as a histogram.K > 1600. Hint: what type of probability distribution would you use for this computation? Investigate the the dramatically useful R function ecdf().K < 1300.K. Hint–use the R quantile() function.# Make a data frame
df = as.data.frame(rbind(zm[[1]], zm[[2]], zm[[3]]))
# Look at the first six rows:
head(df)
## K r sigma tau
## 1 1243.837 0.1982945 0.02899101 1189.7984
## 2 1220.880 0.1987010 0.02724243 1347.4370
## 3 1277.394 0.1869374 0.03609160 767.6931
## 4 1317.190 0.1910025 0.03256109 943.1965
## 5 1272.104 0.2131853 0.03068146 1062.3019
## 6 1267.324 0.2106143 0.03274960 932.3695
# Find the mean of r for the first 1000 iterations
mean(df$r[1: 1000])
## [1] 0.2009405
# Find the mean of r for the last 1000 iterations
nr = length(df$r)
mean(df$r[(nr - 1000): nr])
## [1] 0.2002266
# Make a publication quality plot of the marginal posterior distribution of K as a smooth curve.
# More iterations would produce a smoother cruve) and as a histogram.
plot(density(df$K), main = "", xlim = c(1000, 1500), xlab = "K")
Fig. 2. Posterior density of K.
hist(df$K, main = "", xlim = c(1000, 1500), xlab = "K", breaks = 50, freq = FALSE)
Fig. 2. Posterior density of K.
# Find the probability that the parameter K exceeds 1600
1 - ecdf(df$K)(1600)
## [1] 1e-04
# Find the probability that the parameter 1000 < K < 1300
ecdf(df$K)(1300) - ecdf(df$K)(1000)
## [1] 0.8510333
# Compute the .025 and .975 quantiles of K
quantile(df$K, c(.025, .975))
## 2.5% 97.5%
## 1128.602 1375.907
Summarize the coda output from the logistic model with 4 significant digits. Include Rhat and effective sample size diagnostics (more about these soon). Summarize the coda output for \(r\) alone.
MCMCsummary(zm, round = 4, n.eff = TRUE)
## mean sd 2.5% 50% 97.5% Rhat n.eff
## K 1237.5020 62.7043 1128.6020 1232.8057 1375.9067 1 7016
## r 0.2007 0.0097 0.1813 0.2007 0.2197 1 7235
## sigma 0.0286 0.0030 0.0235 0.0284 0.0353 1 15059
## tau 1259.3378 261.2025 802.0225 1240.1589 1817.1515 1 16330
MCMCsummary(zm, params = "r", round = 4, n.eff = TRUE)
## mean sd 2.5% 50% 97.5% Rhat n.eff
## r 0.2007 0.0097 0.1813 0.2007 0.2197 1 7235
Extract the chains for r and sigma as a data frame using MCMCchains and compute their .025 and .975 quantiles from the extracted object. Display three significant digits. Make a publication quality histogram for the chain for sigma. Indicate the .025 and .975 Bayesian equal-tailed credible interval value with vertical red lines. Overlay the .95 highest posterior density interval with vertical lines in blue. This is a “learn on your own” problem intended to allow you to rehearse the most important goal of this class: being sufficiently confident to figure things out. Load the package HDinterval, enter HDInterval at the console and follow your nose. Also see Hobbs and Hooten Figure 8.2.1.
ex = as.data.frame(MCMCchains(zm, params = c("r", "sigma")))
class(ex)
## [1] "data.frame"
dim(ex)
## [1] 30000 2
signif(apply(ex, 2, quantile, c(.025, .975)), 3)
## r sigma
## 2.5% 0.181 0.0235
## 97.5% 0.220 0.0353
hist(ex$sigma,xlab = expression(sigma), ylab = "Probability density", freq = FALSE, breaks = 100, main = "")
abline(v = quantile(ex$sigma, c(.025, .975)), col = "red", lwd = 2)
abline(v = hdi(ex$sigma, .95), lwd = 2, col = "blue")
For the logistic model:
zm = coda.samples(jm, variable.names = c("K", "r", "sigma", "mu"), n.iter = 10000)
BCI <- MCMCpstr(zm, params = "mu", func = function(x) quantile(x, c(.025, .5, .975)))
HPDI <- MCMCpstr(zm, params = "mu", func = function(x) hdi(x, .95))
plot(data$x, data$y, pch = 19, xlab = "Population size", ylab = "Per-capita grwoth rate (1/year)" )
lines(data$x,BCI$mu[,2], typ = "l")
lines(data$x, BCI$mu[,1], lty = "dashed", col = "red")
lines(data$x, BCI$mu[,3], lty = "dashed", col = "red")
lines(data$x, HPDI$mu[,1], lty = "dashed", col = "blue")
lines(data$x, HPDI$mu[,2], lty = "dashed", col = "blue")
Fig. 4. Median and 95% credible intervals for predicted growth rate and posterior density of K.
The intervals are exactly overlapping in this particular case. Such overlap will not always occur. Equal-tailed intervals based on quantiles will be broader than highest posterior density intervals when marginal posterior distributions are asymmetric. There is a 95% probability that the true value of the mean per-capita population growth rate falls between these lines. Not that this differs from the prediction of a new observation of population growth rate, which would have much broader credible intervals.
Use the MCMCplot function to create caterpillar plots to visualize the posterior densities for \(r\), and \(\sigma\) using the coda object zm from earlier. Then use MCMCplot() to explore different plotting options. There are lots of these options, including ways to make the plots publication quality, show overlap with zero, etc.
MCMCplot(zm, params = "mu")
Rerun the logistic model with n.adapt = 100. Then do the following:
MCMCtrace() and with the Gelman-Rubin, Heidelberger and Welch, and Raftery diagnostics.n.adapt = 100
jm.short = jags.model("LogisticJAGS.R", data = data, inits = inits, n.chains = length(inits), n.adapt = n.adapt)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 3
## Total graph size: 212
##
## Initializing model
n.iter = 500
zm.short = coda.samples(jm.short, variable.names = c("K", "r", "sigma", "tau"), n.iter = n.iter, n.thin = 1)
par(mfrow=c(1,2))
MCMCtrace(zm.short, pdf = FALSE)
gelman.diag(zm.short)
## Potential scale reduction factors:
##
## Point est. Upper C.I.
## K 1.02 1.06
## r 1.01 1.04
## sigma 1.00 1.01
## tau 1.01 1.02
##
## Multivariate psrf
##
## 1.02
heidel.diag(zm.short)
## [[1]]
##
## Stationarity start p-value
## test iteration
## K failed NA 0.00968
## r passed 201 0.13492
## sigma passed 1 0.17658
## tau passed 1 0.17311
##
## Halfwidth Mean Halfwidth
## test
## K <NA> NA NA
## r passed 2.02e-01 2.09e-03
## sigma passed 2.85e-02 3.88e-04
## tau passed 1.28e+03 3.80e+01
##
## [[2]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.558
## r passed 1 0.676
## sigma passed 1 0.547
## tau passed 1 0.587
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 1.58e+01
## r passed 2.00e-01 1.64e-03
## sigma passed 2.88e-02 2.55e-04
## tau passed 1.24e+03 2.14e+01
##
## [[3]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.4864
## r passed 1 0.9910
## sigma passed 101 0.0854
## tau passed 101 0.0795
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 1.17e+01
## r passed 2.01e-01 1.46e-03
## sigma passed 2.85e-02 3.03e-04
## tau passed 1.27e+03 2.71e+01
raftery.diag(zm.short)
## [[1]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
##
## [[2]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
##
## [[3]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
n.update = 500
update(jm.short, n.iter = n.update)
n.iter = 500
zm.short = coda.samples(jm.short, variable.names = c("K", "r", "sigma", "tau"), n.iter = n.iter, n.thin = 1)
MCMCtrace(zm.short, pdf = FALSE)
gelman.diag(zm.short)
## Potential scale reduction factors:
##
## Point est. Upper C.I.
## K 1.01 1.02
## r 1.00 1.01
## sigma 1.00 1.00
## tau 1.00 1.00
##
## Multivariate psrf
##
## 1
heidel.diag(zm.short)
## [[1]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.627
## r passed 1 0.539
## sigma passed 1 0.393
## tau passed 1 0.357
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 2.28e+01
## r passed 2.00e-01 2.94e-03
## sigma passed 2.89e-02 4.48e-04
## tau passed 1.24e+03 3.82e+01
##
## [[2]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.241
## r passed 1 0.395
## sigma passed 1 0.107
## tau passed 1 0.178
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.23e+03 1.57e+01
## r passed 2.01e-01 1.65e-03
## sigma passed 2.87e-02 4.11e-04
## tau passed 1.25e+03 3.49e+01
##
## [[3]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.0824
## r passed 101 0.0941
## sigma passed 101 0.1350
## tau passed 51 0.0599
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 1.41e+01
## r passed 2.00e-01 1.72e-03
## sigma passed 2.85e-02 3.42e-04
## tau passed 1.26e+03 2.68e+01
raftery.diag(zm.short)
## [[1]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
##
## [[2]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
##
## [[3]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
n.update = 10000
update(jm.short, n.iter = n.update)
n.iter = 5000
zm.short = coda.samples(jm.short, variable.names = c("K", "r", "sigma", "tau"), n.iter = n.iter, n.thin = 1)
MCMCtrace(zm.short, pdf = FALSE)
gelman.diag(zm.short)
## Potential scale reduction factors:
##
## Point est. Upper C.I.
## K 1 1.01
## r 1 1.00
## sigma 1 1.00
## tau 1 1.00
##
## Multivariate psrf
##
## 1
heidel.diag(zm.short)
## [[1]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.970
## r passed 1 0.969
## sigma passed 1 0.722
## tau passed 1 0.545
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 5.007758
## r passed 2.01e-01 0.000651
## sigma passed 2.86e-02 0.000118
## tau passed 1.26e+03 9.865859
##
## [[2]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.1000
## r passed 1001 0.0579
## sigma passed 1 0.1397
## tau passed 1 0.0808
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 6.454694
## r passed 2.00e-01 0.000974
## sigma passed 2.86e-02 0.000103
## tau passed 1.26e+03 8.711745
##
## [[3]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.269
## r passed 1 0.930
## sigma passed 1 0.455
## tau passed 1 0.405
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.23e+03 4.310713
## r passed 2.01e-01 0.000570
## sigma passed 2.87e-02 0.000102
## tau passed 1.25e+03 8.254917
raftery.diag(zm.short)
## [[1]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## Burn-in Total Lower bound Dependence
## (M) (N) (Nmin) factor (I)
## K 10 12268 3746 3.27
## r 10 11096 3746 2.96
## sigma 4 5038 3746 1.34
## tau 6 6878 3746 1.84
##
##
## [[2]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## Burn-in Total Lower bound Dependence
## (M) (N) (Nmin) factor (I)
## K 8 8831 3746 2.36
## r 24 23692 3746 6.32
## sigma 5 5391 3746 1.44
## tau 6 6520 3746 1.74
##
##
## [[3]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## Burn-in Total Lower bound Dependence
## (M) (N) (Nmin) factor (I)
## K 10 12122 3746 3.24
## r 4 5299 3746 1.41
## sigma 4 4752 3746 1.27
## tau 5 5673 3746 1.51
Append R code (algorithm 5) to the script you made in exercise 4 to run the JAGS model (algorithm 2) in parallel and estimate the parameters, \(r\), \(K\) \(\sigma\), and \(\tau\). Use the proc.time() function in R to compare the time required for the sequential and parallel JAGS run. If your computer has 3 threads, try running only 2 chains in parallel when doing this exercise. If you have fewer than 3 threads, work with a classmate that has at least 3 threads.
cl <- makeCluster(3)
initsP = function() {list(K = runif(1, 10, 4000), r = runif(1, .01, 2),
sigma = runif(1, 0, 2), .RNG.name = "base::Mersenne-Twister",
.RNG.seed = runif(1, 1, 2000))}
parallel::clusterExport(cl, c("data", "initsP", "n.adapt", "n.update", "n.iter"))
ptm <- proc.time()
out <- clusterEvalQ(cl, {
library(rjags)
jm = jags.model("LogisticJAGS.R", data = data, n.chains = 1,
n.adapt = n.adapt, inits = initsP())
update(jm, n.iter = n.update)
zm = coda.samples(jm, variable.names = c("K", "r", "sigma", "tau"),
n.iter = n.iter, thin = 1)
return(as.mcmc(zm))
})
ParallelTime <- proc.time() - ptm
ParallelTime
## user system elapsed
## 0.004 0.000 0.782
stopCluster(cl)
zmP = mcmc.list(out)
MCMCsummary(zmP)
## mean sd 2.5% 50% 97.5% Rhat
## K 1.237436e+03 6.282117e+01 1.128900e+03 1.232952e+03 1.374828e+03 1
## r 2.007981e-01 9.656524e-03 1.817245e-01 2.008146e-01 2.196020e-01 1
## sigma 2.861924e-02 3.008699e-03 2.349934e-02 2.832982e-02 3.515644e-02 1
## tau 1.260514e+03 2.578836e+02 8.090775e+02 1.245983e+03 1.810876e+03 1
## n.eff
## K 2374
## r 3142
## sigma 9233
## tau 9657
We rerun the model sequentially and use proc.time() again for comparison.
ptm <- proc.time()
jm = jags.model("LogisticJAGS.R", data = data, inits = inits, n.chains = length(inits), n.adapt = n.adapt)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 3
## Total graph size: 212
##
## Initializing model
update(jm, n.iter = n.update)
zm = coda.samples(jm, variable.names = c("K", "r", "sigma"), n.iter = n.iter, n.thin = 1)
SequentialTime <- proc.time() - ptm
SequentialTime
## user system elapsed
## 2.016 0.006 2.024
Looks as if the parallel model runs 2.59 times faster. This factor should increase the more iterations you run (why?).
You might be very tempted, in the name of reproducible science, to set the seed inside each worker to the same value. What happens when you do this? Compare the results from this run to the previous run using MCMCsummary() and MCMC(trace).
cl <- makeCluster(3)
initsP = function() {list(K = runif(1, 10, 4000), r = runif(1, .01, 2),
sigma = runif(1, 0, 2), .RNG.name = "base::Mersenne-Twister",
.RNG.seed = runif(1, 1, 2000))}
parallel::clusterExport(cl, c("data", "initsP", "n.adapt", "n.update", "n.iter"))
out <- clusterEvalQ(cl, {
library(rjags)
set.seed(10)
jm = jags.model("LogisticJAGS.R", data = data, n.chains = 1,
n.adapt = n.adapt, inits = initsP())
update(jm, n.iter = n.update)
zm = coda.samples(jm, variable.names = c("K", "r", "sigma", "tau"),
n.iter = n.iter, thin = 1)
return(as.mcmc(zm))
})
stopCluster(cl)
zmP = mcmc.list(out)
MCMCsummary(zm)
## mean sd 2.5% 50% 97.5% Rhat
## K 1.238707e+03 63.205635682 1130.7480927 1.233152e+03 1.382807e+03 1
## r 2.006284e-01 0.009562590 0.1814912 2.006711e-01 2.192591e-01 1
## sigma 2.867054e-02 0.003024429 0.0235092 2.840384e-02 3.529137e-02 1
## n.eff
## K 1713
## r 2755
## sigma 9365
MCMCsummary(zmP)
## mean sd 2.5% 50% 97.5% Rhat
## K 1.239659e+03 6.227385e+01 1.134518e+03 1.233489e+03 1378.0986189 NaN
## r 2.005185e-01 9.584991e-03 1.817614e-01 2.004705e-01 0.2194044 NaN
## sigma 2.859709e-02 2.896285e-03 2.358206e-02 2.843383e-02 0.0348032 NaN
## tau 1.260124e+03 2.518494e+02 8.255847e+02 1.236885e+03 1798.1932583 NaN
## n.eff
## K 3953
## r 4096
## sigma 9014
## tau 10309
MCMCtrace(zm, pdf = FALSE)
MCMCtrace(zmP, pdf = FALSE)
Repeat this exercise with fixed initial values using algorithm 6.
cl <- makeCluster(3)
myWorkers <- NA
for(i in 1:3) myWorkers[i] <- word(capture.output(cl[[i]]),-1)
initsP = list(
list(K = 1500, r = .2, sigma = 1, RNG.seed = 1,
.RNG.name = "base::Mersenne-Twister"),
list(K = 1000, r = .15, sigma = .1, RNG.seed = 23,
.RNG.name = "base::Mersenne-Twister"),
list(K = 900, r = .3, sigma = .01, .RNG.seed = 255,
.RNG.name = "base::Mersenne-Twister"))
parallel::clusterExport(cl, c("myWorkers", "data", "initsP", "n.adapt", "n.update", "n.iter"))
out <- clusterEvalQ(cl, {
library(rjags)
jm = jags.model("LogisticJAGS.R", data = data, n.chains = 1,
n.adapt = n.adapt, inits = initsP[[which(myWorkers==Sys.getpid())]])
update(jm, n.iter = n.update)
zm = coda.samples(jm, variable.names = c("K", "r", "sigma", "tau"),
n.iter = n.iter, thin = 1)
return(as.mcmc(zm))
})
stopCluster(cl)
zmP = mcmc.list(out)
MCMCsummary(zmP)
## mean sd 2.5% 50% 97.5% Rhat
## K 1.236654e+03 6.207641e+01 1.129370e+03 1.232326e+03 1375.7625637 1.01
## r 2.009182e-01 9.614570e-03 1.817771e-01 2.009465e-01 0.2197338 1.00
## sigma 2.863496e-02 3.043586e-03 2.351144e-02 2.840583e-02 0.0353045 1.00
## tau 1.260038e+03 2.606176e+02 8.023057e+02 1.239324e+03 1809.0121221 1.00
## n.eff
## K 1715
## r 2558
## sigma 7247
## tau 7800
MCMCtrace(zmP, pdf = FALSE)